
#include "Self_Define_Functions.h"

using namespace std;

//compare two struct variables
bool Compare_NewType( NewType data1, NewType data2 )
{
    return data1.data > data2.data;
}

//compute the fitness of an individual
double Fitness ( double *particle, int &Fes, CEC2013 *pFunc )
{
    Fes++;
    return pFunc->evaluate( particle );
}

//compute the fitness of a population
void Popupation_Fitness( double **population, int population_size, int &Fes, double *result, CEC2013 *pFunc )
{
    int i;
    for ( i = 0; i < population_size; ++i )
    {
        result[i] = pFunc->evaluate( population[i] );
    }

    Fes += population_size;
}

//compute the Euclidean distance between two vectors
double Distance( double *vector1, double *vector2, int dim )
{
    double sum = 0;
    for( int i = 0; i < dim; ++i )
    {
        sum += pow( vector1[i] - vector2[i], 2);
    }

    return sqrt( sum );
}

//compute the Euclidean distance between two vectors
double Distance( double *vector1, vector<double> vector2, int dim )
{
    assert( vector2.size() == dim );

    double sum = 0;
    for( int i = 0; i < dim; ++i )
    {
        sum += pow( vector1[i] - vector2[i], 2);
    }

    return sqrt( sum );
}

//obtain the candidate global optima of a population
void Get_Seeds( double **population, double *population_results, int population_size, int dim, vector<double> &seed_fitness, double radius )
{
    int i,j;
    bool found;
    double dist;

    vector< vector<double> > seeds;

    //sort population
    NewType *temp = new NewType [population_size];
    for( i = 0; i < population_size; ++i )
    {
        temp[i].data = population_results[i];
        temp[i].id = i;
    }


    sort( temp, temp+population_size, Compare_NewType );

    seed_fitness.clear();

    for( i = 0; i < population_size; ++i )
    {
        found = false;

        for( j = 0; j < seeds.size(); ++j )
        {
            dist = Distance( population[temp[i].id], seeds[j], dim );
            if( dist <= radius )
            {
                found = true;
                break;
            }
        }

        if( !found )
        {
            vector<double> temp_seed( population[temp[i].id], population[temp[i].id] + dim );
            seeds.push_back(temp_seed);
            seed_fitness.push_back( population_results[temp[i].id] );
        }
    }


    delete []temp;
    seeds.clear();
}

//obtain the number of global optima
int How_Many_Global_Optima( vector<double> seed_fitness, double epsilon, CEC2013 *pFunc )
{

    int  counter = 0;;
	for ( unsigned int  i = 0; i < seed_fitness.size(); ++i )
	{
		/* |F_seed - F_goptimum| <= accuracy */
		if ( fabs(seed_fitness[i] - pFunc->get_fitness_goptima()) <= epsilon)
		{
			++counter;
		}
		/* save time */
		if (counter == pFunc->get_no_goptima())
            break;
	}
	return counter;

}

//obtain the number of found global optima in a population
int Compute_Global_Optima_Found( double **population, double *population_results, int population_size, int dim, double epsilon, CEC2013 *pFunc )
{
    vector<double> seed_fitness;

    Get_Seeds( population, population_results,  population_size, dim, seed_fitness, pFunc->get_rho() );

    int counter = How_Many_Global_Optima( seed_fitness,  epsilon, pFunc );

    seed_fitness.clear();

    return counter;

}

//find the closest individual of the child
int Close_Particle( double *child, double **population, int population_size, int dim )
{
    double min_distance = 1e100;

    int i;

    int particle_index = -1;

    double temp_distance;

    for( i = 0; i < population_size; ++i )
    {
        temp_distance = Distance( child, population[i], dim );
        if( min_distance > temp_distance )
        {
            min_distance = temp_distance;
            particle_index = i;
        }
    }

    return particle_index;
}

//find the closest individual based on the computed distance
int Close_Particle( double *dist,  int population_size )
{
    double min_distance = dist[0];

    int i;

    int particle_index = 0;


    for( i = 1; i < population_size; ++i )
    {
        if( min_distance > dist[i] )
        {
            min_distance = dist[i];
            particle_index = i;
        }
    }

    return particle_index;
}


void Dist_Individual_to_All( double *dist, double *individual, double **population, int population_size, int dim )
{
    int i;
    for( i = 0; i < population_size; ++i )
    {
        dist[i] = Distance( individual, population[i], dim );
    }

}

//compute the distance between an individual and each individual in the  population
void Dist_All_to_All( double **population, double **dist, int population_size, int dim )
{
    int i,j;
    for( i = 0; i < population_size; ++i )
    {
        dist[i][i] = 0;
    }

    for( i = 0; i < population_size; ++i )
    {
        for( j = i+1; j < population_size; ++j )
        {
            dist[i][j] = dist[j][i] = Distance( population[i], population[j], dim );
        }
    }
}

//update the distance matrix based on a distance vector
void Update_Dist( double **dist, double *temp_dist, int closet_individual, int population_size )
{
    int i;

    memcpy( dist[closet_individual], temp_dist, sizeof(double)*population_size );

    for( i = 0; i < population_size; ++i )
        dist[i][closet_individual] = temp_dist[i];
}

//sort the population based on the computed distance
void Sort_Population_Dist( double *dist, int *sorted_index, int population_size )
{
    int i;

    NewType *temp = new NewType [population_size];

    for( i = 0; i < population_size; ++i )
    {
        temp[i].data = dist[i];
        temp[i].id = i;
    }


    sort( temp, temp+population_size, Compare_NewType );

    for( i = 0; i < population_size; ++i )
        sorted_index[i] = temp[population_size-1-i].id;

    delete []temp;
}

//sort the population based on the computed fitness
void Sort_Population_Fitness( double *population_results, int *sorted_index, int population_size )
{
    int i;

    NewType *temp = new NewType [population_size];

    for( i = 0; i < population_size; ++i )
    {
        temp[i].data = population_results[i];
        temp[i].id = i;
    }

    sort( temp, temp+population_size, Compare_NewType );

    for( i = 0; i < population_size; ++i )
        sorted_index[i] = temp[i].id;

    delete []temp;
}

//obtain the ranking of each individual based on the computed fitness
void Ranking( double *population_results, int *ranking, int population_size, double &fitness_gap )
{
    int i;

    NewType *temp = new NewType [population_size];

    for( i = 0; i < population_size; ++i )
    {
        temp[i].data = population_results[i];
        temp[i].id = i;
    }

    sort( temp, temp+population_size, Compare_NewType );

    for( i = 0; i < population_size; ++i )
        ranking[temp[i].id] = i;

    fitness_gap = temp[0].data - temp[population_size-1].data;

    delete []temp;
}

//obtain the M nearest individuals of an individual based on the computed distance
void Get_Nearest_M_Individuals( double *dist, bool *flag, int *Nearest_M, int population_size, int n_nearest )
{
    int i,j;

    NewType *temp = new NewType [population_size];
    for( i = 0; i < population_size; ++i )
    {
        temp[i].data = dist[i];
        temp[i].id = i;
    }

    sort( temp, temp+population_size, Compare_NewType );

    j = population_size-1;
    for( i = 1; i < n_nearest; ++i )
    {
        while( j >=0 )
        {
            if( !flag[temp[j].id] )
            {
                Nearest_M[i] = temp[j].id;
                flag[temp[j].id] = true;
                --j;
                break;
            }
            else
            {
                --j;
            }

        }

    }

    delete []temp;

}

//divide the population into niches (crowding)
void Species( double *reference_point, double **population, int **species, double **dist,
              int species_num, int N_nearest, int last_species_size, int population_size, int dim )
{
    int i,j;
    int start_search = 0;
    int final_size;

    bool *flag = new bool[population_size];
    int *sorted_index = new int [population_size];

    for( i = 0; i < population_size; ++i )
        flag[i] = false;

    double *ref_dist = new double[population_size];

    //compute the Euclidean distance between the reference point and the population
    Dist_Individual_to_All( ref_dist, reference_point, population, population_size, dim );

    //sort the distance
    Sort_Population_Dist( ref_dist, sorted_index, population_size );

    start_search = 0;

    for( i = 0; i < species_num-1; ++i )
    {

        species[i][0] = sorted_index[start_search];//find seed
        flag[sorted_index[start_search]] = true;

        //get the M nearest individuals
        Get_Nearest_M_Individuals( dist[species[i][0]], flag, species[i], population_size, N_nearest );

        //update start_search
        while( flag[sorted_index[start_search]] )
            ++start_search;
    }

    if( last_species_size == 0 )
    {
        final_size = N_nearest;
    }
    else
    {
        final_size = last_species_size;
    }

    for( j = 0; j < final_size; ++j )
    {
         while( start_search < population_size )
        {
            if( !flag[sorted_index[start_search]] )
            {
                species[i][j] = sorted_index[start_search];
                 flag[sorted_index[start_search]] = true;
                 ++start_search;
                 break;

            }
            else
            {
                ++start_search;
            }
        }

    }

    delete []ref_dist;
    delete []sorted_index;
    delete []flag;
}

//obtain the seed of each niche
void Get_Seed_Location( int *seed_location, double *population_results, int **species, int species_num, int speices_size, int last_species_size )
{
    int i,j;
    double best;
    int final_size;

    //find the seed of the 1st to species_num-1 th niche
    for( i = 0; i < species_num-1; ++i )
    {
        best = population_results[species[i][0]];
        seed_location[i] = species[i][0];

        for( j = 1; j < speices_size; ++j )
        {
            if( best < population_results[species[i][j]] )
            {
                best = population_results[species[i][j]];
                seed_location[i] = species[i][j];
            }
        }
    }

    if( last_species_size == 0 )
    {
        final_size = speices_size;
    }
    else
    {
        final_size = last_species_size;
    }

    //find the seed of the last niche
    best = population_results[species[i][0]];
    seed_location[i] = species[i][0];

    for( j = 1; j < final_size; ++j )
    {
        if( best < population_results[species[i][j]] )
        {
            best = population_results[species[i][j]];
            seed_location[i] = species[i][j];
        }
    }
}

//compute the probability of each niche to conduct local search
void Seed_Local_Prob( double *seed_local_prob, int *seed_location, double *population_results, int species_num )
{
    int i;

    double min_value;
    double max_value;
    double normalize_factor;

    min_value = max_value = population_results[seed_location[0]];

    //find the max and min fitness of the seed set
    for( i = 1; i < species_num; ++i )
    {
        if( min_value > population_results[seed_location[i]] )
            min_value = population_results[seed_location[i]];
        if( max_value < population_results[seed_location[i]] )
            max_value = population_results[seed_location[i]];
    }

    if( min_value <= 0 )
    {
        normalize_factor = max_value + fabs(min_value)+1;
        for( i = 0; i < species_num; ++i )
        {
            seed_local_prob[i] = ( population_results[seed_location[i]] + fabs(min_value) +1 ) / normalize_factor;
        }
    }
    else
    {
        for( i = 0; i < species_num; ++i )
        {
            seed_local_prob[i] =  population_results[seed_location[i]] / max_value;
        }
    }
}

//local search conducted on the seed of each niche
void Local_Species_Evolve(  double **population, double *population_results, double *LBound, double *UBound, int seed_location,
                            int dim,  double local_std_value , int &Fes, CEC2013 *pFun, int sample_num )
{
    int j,k;
    boost::mt19937 generator(time(0)*rand());

    double *temp_particle = new double[dim];
    double temp_fitness ;


        for( k = 0; k < sample_num; ++k )
        {
            for( j = 0; j < dim; ++j )
            {

                boost::normal_distribution<> norm_generator( population[seed_location][j], local_std_value );
                boost::variate_generator< boost::mt19937&, boost::normal_distribution<> > norm_generation( generator, norm_generator );


                temp_particle[j] = norm_generation();
                if( temp_particle[j] < LBound[j] )
                    temp_particle[j] = LBound[j];

                if( temp_particle[j] > UBound[j] )
                    temp_particle[j] = UBound[j];
            }

            temp_fitness = pFun->evaluate(temp_particle);
            Fes++;
            if( temp_fitness > population_results[seed_location] )
            {
                population_results[seed_location] = temp_fitness;
                memcpy( population[seed_location], temp_particle, sizeof(double)*dim );
            }

        }

    delete []temp_particle;
}

//local search conducted on the seed of each niche
void Local_Evolve(  double **population, double *population_results, double *LBound, double *UBound, double *seed_local_prob, int *seed_location,
                    int dim, int species_num,  double local_std_value, int &Fes, CEC2013 *pFun, int sample_num )
{
    int i;

    double temp_prob;

    boost::mt19937 generator(time(0)*rand());
    boost::uniform_real<> uniform_real_generate_r( 0, 1 );
    boost::variate_generator< boost::mt19937&, boost::uniform_real<> > random_real_num_r( generator, uniform_real_generate_r );

    for( i = 0; i < species_num; ++i )
    {
        temp_prob = random_real_num_r();
        if( temp_prob <= seed_local_prob[i] )
            Local_Species_Evolve(  population, population_results, LBound, UBound, seed_location[i], dim, local_std_value, Fes, pFun, sample_num );
    }
}

//compute the weight of each ant in each species
void Weight_in_Species( double *weight, int *species, double*population_results, int species_size, double population_fitness_gap )
{
    int i;

    double weight_parameter;

    double species_fitness_gap = 0;

    int *ranking = new int[species_size];

    double*species_results = new double[species_size];

    for( i =0; i < species_size; ++i )
    {
        species_results[i] = population_results[species[i]];
    }

    Ranking( species_results, ranking, species_size, species_fitness_gap );

    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    if( population_fitness_gap != 0 )
        weight_parameter = 0.1 + 0.3*exp(-1*(species_fitness_gap/population_fitness_gap));
    else
        weight_parameter = 0.1 + 0.3;
    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////


    for( i = 0; i < species_size; ++i )
    {
        weight[i] = 1.0/( weight_parameter * species_size * sqrt( 2*M_PI) )* exp( -1*(ranking[i]*ranking[i])/( 2*weight_parameter*weight_parameter*species_size*species_size ) );
    }

    delete []ranking;
    delete []species_results;
}

//compute the selection probability of each ant in each species
void Probability_in_Species( double *probability,  double *weight, int species_size )
{
    int i;
    double sum_weight = 0;
    for( i = 0; i < species_size; ++i )
    {
        sum_weight += weight[i];
    }

    for( i = 0; i < species_size; ++i )
    {
        probability[i] = weight[i] / sum_weight;
    }

    for( i = 1; i < species_size; ++i )
    {
        probability[i] += probability[i-1];
    }
}


double Cal_Variance( int *species, double **population, int dim_index, double dim_value, double epsilon_norm, int species_size )
{
    int i;
    double variance = 0;
    for( i = 0; i < species_size; ++i )
    {
        variance += abs( population[species[i]][dim_index] - dim_value );
    }

    variance = epsilon_norm * variance /(species_size-1);

    if( variance == 0 )
        variance = 1e-4;

    return variance;
}


void Cal_Variance( double *variance, int *species, double **population, double *selected_individual,  double epsilon_norm, int species_size, int dim )
{
    int i;

    if( species_size == 1 )
    {
        for( i = 0; i < dim; ++i )
        {
            variance[i] = 1e-4;
        }
    }
    else
    {

        for( i = 0; i < dim; ++i )
        {
            variance[i] = Cal_Variance( species, population, i, selected_individual[i], epsilon_norm, species_size );
        }
    }
}


//randomly select an individual
int Selected_Individual( double *probability, int species_size )
{
    int i;
    boost::mt19937 generator(time(0)*rand());
    boost::uniform_real<> uniform_real_generate_r( 0, 1 );
    boost::variate_generator< boost::mt19937&, boost::uniform_real<> > random_real_num_r( generator, uniform_real_generate_r );

    double random_value = random_real_num_r();

    if( random_value >=0 && random_value <= probability[0] )
        return 0;

    for( i = 1; i < species_size; ++i )
    {
        if( random_value > probability[i-1] && random_value < probability[i] )
            break;
    }

    return i;
}


//evolve each species
void Species_Evolve_At_Individual3( int *species, double **population, double*population_results, double **child, double *LBound, double *UBound, int seed_location,
                                    int dim, int species_size, double population_fitness_gap )
{
    int i,j;
    int selected_index;
    double *std_value = new double[dim];

    double *weight = new double[species_size];
    double *probability = new double[species_size];
    double *mean = new double[dim];
    double F;

    double epsilon_norm;


    boost::mt19937 generator(time(0)*rand());
    boost::uniform_real<> uniform_real_generator2( 0, 1 );
    boost::variate_generator< boost::mt19937&, boost::uniform_real<> > random_real_param( generator, uniform_real_generator2 );


    boost::uniform_real<> uniform_real_generator( 0, 1 );
    boost::variate_generator< boost::mt19937&, boost::uniform_real<> > random_real_epsilon_norm( generator, uniform_real_generator );

    //compute the weight of each solution in the species
    Weight_in_Species( weight, species, population_results, species_size, population_fitness_gap );

    //compute the selection probability of each solution in the species
    Probability_in_Species( probability, weight, species_size );

    for( j = 0; j < species_size; ++j )
    {
        //randomly select a solution
        selected_index = Selected_Individual( probability, species_size );
        if( random_real_param() < 0.5 )
            memcpy( mean, population[species[selected_index]], sizeof(double)*dim );//use the selected individual as the mean
        else//use DE to mutate the mean
        {
            F = random_real_param();
            for( i = 0; i < dim; ++i )
            {
                mean[i] = population[species[selected_index]][i] + F*(population[seed_location][i] - population[species[selected_index]][i] );
            }
        }
        epsilon_norm  = random_real_epsilon_norm();

        //compute the variance
        Cal_Variance( std_value, species, population, mean, epsilon_norm, species_size, dim );


        for( i = 0; i < dim; ++i )//randomly generate an offspring
        {
            boost::normal_distribution<> norm_generator( mean[i], std_value[i] );
            boost::variate_generator< boost::mt19937&, boost::normal_distribution<> > norm_generation( generator, norm_generator );
            child[species[j]][i] = norm_generation();
            if( child[species[j]][i] < LBound[i] )
                child[species[j]][i] = LBound[i];

            if( child[species[j]][i] > UBound[i] )
                child[species[j]][i] = UBound[i];
        }

    }

    delete []std_value;
    delete []weight;
    delete []probability;
    delete []mean;
}



//evolve the population
void Evolve_At_Individual3( int **species, double **population, double *population_results, double **child, double *LBound, double *UBound, int *seed_location,
                            int population_size, int dim, int species_num, int n_nearest, int last_species_size )
{
    int i;

    double population_fitness_gap;

    population_fitness_gap = Get_Population_Fitness_Gap( population_results, population_size );

    for( i = 0; i < species_num-1; ++i )//evolve the 1st to species_num-1 th species
    {
        Species_Evolve_At_Individual3( species[i], population, population_results, child, LBound, UBound,seed_location[i], dim, n_nearest, population_fitness_gap );
    }

    if( last_species_size == 0 )//evolve the last species
        Species_Evolve_At_Individual3( species[i], population, population_results, child, LBound, UBound,seed_location[i], dim, n_nearest, population_fitness_gap );
    else
        Species_Evolve_At_Individual3( species[i], population, population_results, child, LBound, UBound,seed_location[i], dim, last_species_size, population_fitness_gap );

}

//get f_max-f_min namely the best fitness - the worst fitness
double Get_Population_Fitness_Gap( double *population_results, int population_size )
{
    double fitness_min, fitness_max;

    int i;

    fitness_max = fitness_min = population_results[0];

    for( i = 1; i < population_size; ++i )
    {
        if( fitness_max < population_results[i] )
            fitness_max = population_results[i];

        if(fitness_min > population_results[i] )
            fitness_min = population_results[i];

    }

    return (fitness_max-fitness_min);

}

//select individuals for next generation
void Selection( double **population, double **child, double **dist, double *population_results, double *child_results,
                int population_size, int dim )
{
    int i;

    int closest_index;

    double *temp_dist = new double[population_size];

    for( i =0 ; i < population_size; ++i )
    {
        Dist_Individual_to_All( temp_dist, child[i], population, population_size, dim );

        closest_index = Close_Particle( temp_dist, population_size );

        if( child_results[i] > population_results[closest_index] )
        {
            memcpy( population[closest_index], child[i], sizeof(double)*dim );
            population_results[closest_index] = child_results[i];
            temp_dist[closest_index] = 0;
            Update_Dist( dist, temp_dist, closest_index, population_size );
        }
    }

    delete []temp_dist;
}

//randomly generate a reference point
void Generate_Reference_Point( double *reference_point, double *LBound, double *UBound, int dim )
{
    int j;
    for( j =0; j < dim; ++j )
    {
        boost::mt19937 generator(time(0)*rand());
        boost::uniform_real<> uniform_real_generate_x( LBound[j], UBound[j] );
        boost::variate_generator< boost::mt19937&, boost::uniform_real<> > random_real_num_x( generator, uniform_real_generate_x );

        reference_point[j] = random_real_num_x();
    }
}

//update the distance  between individuals
void Update_Dist_Seed( double **dist, int *seed_location, double **population, int population_size, int seed_num, int dim )
{
    double *temp_individual = new double[dim];
    double *temp_dist = new double[population_size];

    int i;

    for( i = 0; i < seed_num; ++i )
    {
        memcpy( temp_individual, population[seed_location[i]], sizeof(double)*dim );
        Dist_Individual_to_All( temp_dist, temp_individual, population, population_size, dim );
        temp_dist[seed_location[i]] = 0;
        Update_Dist( dist, temp_dist, seed_location[i], population_size );
    }

    delete []temp_individual;
    delete []temp_dist;
}

































